# take a look at the data
glimpse(med_roads)
head(med_roads)
# get rid of everything except the gis field
med_roads <- med_roads %>% dplyr::select()
# crop the roads
med_roads <- med_roads %>%
st_crop(bbox_cities)
glimpse(med_roads)
# add Via Reggia
crs = st_crs(4326)
s1 <- rbind(c(24.728611,49.124722),c(14.4378005, 50.0755381))
ls <- st_linestring(s1)
b = st_sf(geometry = st_sfc(ls), crs = crs)
med_roads <- med_roads %>%
rbind(b)
# plot the roads
plot(med_roads[1], col="brown", reset=F, add=T)
# Bring in rivers shape file: round 1
rivers <- st_read("Rivers/MajorRivers.shp") %>%
st_transform(4326)
glimpse(rivers)
rivers <- rivers %>% dplyr::select(NAME)
rivers <- rivers %>%
st_crop(bbox_cities)
glimpse(rivers)
plot(rivers[2], col="blue", reset=F, add=T)
# Bring in rivers shape file: round 2
rivers2 <- st_read("Rivers/Rivers_conic.shp") %>%
st_transform(4326)
glimpse(rivers2)
rivers2 <- rivers2 %>% dplyr::select(name)
rivers2 <- rivers2 %>%
st_crop(bbox_cities)
glimpse(rivers2)
target <- c("Dniester", "Volga", "Oka")
rivers2 <- rivers2 %>% filter(name %in% target)
plot(rivers2[2], col="blue", reset=F, add=T)
# Bring in the modern country borders
modern_countries <- st_read("ModernCountries/Longitude_Graticules_and_World_Countries_Boundaries-shp/99bfd9e7-bb42-4728-87b5-07f8c8ac631c2020328-1-1vef4ev.lu5nk.shp") %>%
st_transform(4326)
# Clip the countries using bounding box based on the cities data
modern_countries <- modern_countries %>%
st_crop(bbox_cities)
glimpse(modern_countries)
head(modern_countries)
modern_countries <- modern_countries %>% dplyr::select(CNTRY_NAME)
glimpse(modern_countries)
plot(modern_countries[2], col="NA", reset=F, add=T)
# Bring in seas shape file
seas <- st_read("Seas/ne_10m_ocean.shp") %>%
st_transform(4326)
glimpse(seas)
# there is a ring intersection error in the seas .shp file. So need to fix the geometry.
sf::sf_use_s2(FALSE) # this does work
seas_crop <- st_crop(seas, bbox_cities)
sf::sf_use_s2(TRUE) # putting settings back to their default
# make figure A1
state_map <- tm_shape(modern_countries) + tm_borders(col = "grey")
rivers_map <- tm_shape(rivers) +
tm_lines(col="blue", lwd=3) +
tm_add_legend("line", col = "blue", title = "Major Rivers", lwd=3)
rivers_map2 <- tm_shape(rivers2) +
tm_lines(col="blue", lwd=3)
med_map <- tm_shape(med_roads) +
tm_lines(col="black") +
tm_add_legend("line", col = "black", title = "Trade Routes")
seas_map <- tm_shape(seas) +
tm_polygons(col="grey")
cities_map <- tm_shape(cities) +
tm_text("city_jjk", size=0.75,ymod=1,xmod=.75, fontface="bold") +
tm_dots(col="red",size=.40) +
tm_add_legend("symbol", col = "red", title = "Trade Cities")
full_map <- rivers_map + rivers_map2 + med_map + cities_map + seas_map
full_map <- full_map + tm_layout(legend.position = c(.77,.70),legend.frame=T)
full_map
# Code to create Figure A1 in "Pandemics and Cities: Evidence from the Black Death and the Long-Run"
# Author: Noel Johnson
# Date Created: 5-14-20
# Last updated: 9-5-23
library(sf)
library(raster)
library(tidyverse)
library(tmap)
library(haven)
# Need the following features in the map:
# (1) Modern Country Borders
# (2) Cities: ASTRAHAN, FEODOSIJA, MESSINA, GENOVA, IZMAIL, WIEN, MOSKVA, NOVGOROD, VISBY, Prague, Leipzig, Halych
# (3) Rivers: Danube, Dniester, the river from Astrakhan to Moscow, moscow to Novgord
# (4) Medieval Trade Routes: the medieval land route (Via Regia),
# (5) seas
# coordinates of the cities not in Bosker (from wikipedia)
# c("ASTRAHAN",51.166667, 71.433333)
# c("FEODOSIJA",45.048889, 35.379167)
# c("VINICA",45.351667, 28.836389)
# c("MOSKVA",55.755833, 37.617222)
# c("NOVGOROD",58.55, 31.266667)
# c("VISBY",57.634722, 18.299167)
# set working directory
setwd("/Users/noeljohnson_laptop/Dropbox/Research/Plague City Growth/JUE RR1 Nov 2022/figures_replication/Figure_A1")
# bring in the cities
cities_raw <- read_dta("All_Cities.dta")
# keep the variables you want
cities_raw <- cities_raw %>% select(city_jjk, longitude, latitude)
# keep the cities you want
target <- c("WIEN", "PRAHA", "LEIPZIG")
cities_raw <- cities_raw %>% filter(city_jjk %in% target)
# add the cities not contained in Bosker that you acquired coordinates for form Bosker
cities_raw <- cities_raw %>%
rbind(c("ASTRAHAN (1345)",48.05,46.35)) %>%
rbind(c("KAFFA (1346)",35.379167,45.048889)) %>%
rbind(c("VINICA",28.836389,45.351667)) %>%
rbind(c("MOSKVA",37.617222,55.755833)) %>%
rbind(c("NOVGOROD",31.266667,58.55)) %>%
rbind(c("VISBY",18.299167,57.634722)) %>%
rbind(c("HALYCH",24.728611,49.124722)) %>%
rbind(c("GENOVA",8.946256,44.4056499)) %>%
rbind(c("MESSINA (1347)",15.5540152,38.1938137))
# clean up the cities data
glimpse(cities_raw)
cities_raw$longitude <- as.numeric(as.character(cities_raw$longitude))
cities_raw$latitude <- as.numeric(as.character(cities_raw$latitude))
glimpse(cities_raw)
# spatialize the cities data
cities <- st_as_sf(cities_raw, coords=c("longitude", "latitude"), crs=4326)
# take a look at the data
cities
glimpse(cities)
head(cities)
plot(cities[2], col="red", pch=16, reset=F)
# Make a bounding box of the cities so you can create one a little bigger
bbox_cities <- st_bbox(cities, crs=4326)
bbox_cities
xrange <- bbox_cities$xmax - bbox_cities$xmin # range of x values
yrange <- bbox_cities$ymax - bbox_cities$ymin # range of y values
bbox_cities[1] <- bbox_cities[1] - (0.20 * xrange) # xmin - left
bbox_cities[3] <- bbox_cities[3] + (0.20 * xrange) # xmax - right
bbox_cities[2] <- bbox_cities[2] - (0.20 * yrange) # ymin - bottom
bbox_cities[4] <- bbox_cities[4] + (0.20 * yrange) # ymax - top
# crop the cities using your bounding box
cities <- cities %>%
st_crop(bbox_cities)
# Bring in medieval roads shape file
med_roads <- st_read("Europe_Medieval_Trade_Routes/MedRdsProj.shp") %>%
st_transform(4326)
# take a look at the data
glimpse(med_roads)
head(med_roads)
# get rid of everything except the gis field
med_roads <- med_roads %>% dplyr::select()
# crop the roads
med_roads <- med_roads %>%
st_crop(bbox_cities)
glimpse(med_roads)
# add Via Reggia
crs = st_crs(4326)
s1 <- rbind(c(24.728611,49.124722),c(14.4378005, 50.0755381))
ls <- st_linestring(s1)
b = st_sf(geometry = st_sfc(ls), crs = crs)
med_roads <- med_roads %>%
rbind(b)
# plot the roads
plot(med_roads[1], col="brown", reset=F, add=T)
# Bring in rivers shape file: round 1
rivers <- st_read("Rivers/MajorRivers.shp") %>%
st_transform(4326)
glimpse(rivers)
rivers <- rivers %>% dplyr::select(NAME)
rivers <- rivers %>%
st_crop(bbox_cities)
glimpse(rivers)
plot(rivers[2], col="blue", reset=F, add=T)
# Bring in rivers shape file: round 2
rivers2 <- st_read("Rivers/Rivers_conic.shp") %>%
st_transform(4326)
glimpse(rivers2)
rivers2 <- rivers2 %>% dplyr::select(name)
rivers2 <- rivers2 %>%
st_crop(bbox_cities)
glimpse(rivers2)
target <- c("Dniester", "Volga", "Oka")
rivers2 <- rivers2 %>% filter(name %in% target)
plot(rivers2[2], col="blue", reset=F, add=T)
# Bring in the modern country borders
modern_countries <- st_read("ModernCountries/Longitude_Graticules_and_World_Countries_Boundaries-shp/99bfd9e7-bb42-4728-87b5-07f8c8ac631c2020328-1-1vef4ev.lu5nk.shp") %>%
st_transform(4326)
# Clip the countries using bounding box based on the cities data
modern_countries <- modern_countries %>%
st_crop(bbox_cities)
glimpse(modern_countries)
head(modern_countries)
modern_countries <- modern_countries %>% dplyr::select(CNTRY_NAME)
glimpse(modern_countries)
plot(modern_countries[2], col="NA", reset=F, add=T)
# Bring in seas shape file
seas <- st_read("Seas/ne_10m_ocean.shp") %>%
st_transform(4326)
glimpse(seas)
# there is a ring intersection error in the seas .shp file. So need to fix the geometry.
sf::sf_use_s2(FALSE) # this does work
seas_crop <- st_crop(seas, bbox_cities)
sf::sf_use_s2(TRUE) # putting settings back to their default
# make figure A1
state_map <- tm_shape(modern_countries) + tm_borders(col = "grey")
rivers_map <- tm_shape(rivers) +
tm_lines(col="blue", lwd=3) +
tm_add_legend("line", col = "blue", title = "Major Rivers", lwd=3)
rivers_map2 <- tm_shape(rivers2) +
tm_lines(col="blue", lwd=3)
med_map <- tm_shape(med_roads) +
tm_lines(col="black") +
tm_add_legend("line", col = "black", title = "Trade Routes")
seas_map <- tm_shape(seas_crop) +
tm_polygons(col="grey")
cities_map <- tm_shape(cities) +
tm_text("city_jjk", size=0.75,ymod=1,xmod=.75, fontface="bold") +
tm_dots(col="red",size=.40) +
tm_add_legend("symbol", col = "red", title = "Trade Cities")
full_map <- rivers_map + rivers_map2 + med_map + cities_map + seas_map
full_map <- full_map + tm_layout(legend.position = c(.77,.70),legend.frame=T)
full_map
# save the map (note: the orange arrows were inputted by hand)
tmap_save(full_map, filename="Trade_Figure1.pdf",
height=8.5, width=11, units="in", dpi=300)
# Code to create Figure A5 in "Pandemics and Cities: Evidence from the Black
# Death and the Long-Run"
# Author: Noel Johnson
# Date Created: 12-10-21
# Last updated: 9-8-23
library(tidyverse)
library(sf)
library(tmap)
library(haven)
# set working directory
setwd("/Users/noeljohnson_laptop/Dropbox/Research/Plague City Growth/JUE RR1 Nov 2022/figures_replication/Figure_A5")
# read in, project, and plot the 1300 borders map from Nussli
map_1300 <- st_read("gis_1300/data/1300/latin1/sovereign_dependent_states.shp")
map_1300 <- st_transform(map_1300, crs=4326)
map_1300
glimpse(map_1300)
plot(map_1300[12], axes=T, reset=T)
# read in, project, and plot the cities
cities_165 <- read_csv("cities_for_map165.csv")
# spatialize the city data. turn into sf points, specify coords and crs
cities_165_sf <- st_as_sf(cities_165,
coords = c("longitude", "latitude"),
crs = 4326)
plot(cities_165_sf[4], axes=T, reset=T)
# plot the borders and cities together
plot(map_1300[12], reset=F)
plot(cities_165_sf[4], add=T)
# need to crop the map and make it look prettier
# Make a bounding box of the cities so you can create one a little bigger
bbox_cities <- st_bbox(cities_165_sf, crs=4326)
bbox_cities
xrange <- bbox_cities$xmax - bbox_cities$xmin # range of x values
yrange <- bbox_cities$ymax - bbox_cities$ymin # range of y values
bbox_cities[1] <- bbox_cities[1] - (0.1 * xrange) # xmin - left
bbox_cities[3] <- bbox_cities[3] + (0.1 * xrange) # xmax - right
bbox_cities[2] <- bbox_cities[2] - (0.1 * yrange) # ymin - bottom
bbox_cities[4] <- bbox_cities[4] + (0.1 * yrange) # ymax - top
bbox_cities_sf <- bbox_cities %>%  # take the bounding box ...
st_as_sfc() # ... and make it a sf polygon
bbox_cities_sf
# plot(bbox_cities_sf)
# plot the borders, cities, and bounding box together
plot(map_1300[12], reset=F)
plot(cities_165_sf[4], add=T)
plot(bbox_cities_sf, add=T)
# clip the map and cities using the bounding box
cities_165_crop <- st_crop(cities_165_sf, bbox_cities_sf)
# have to fix the geometry of map_1300 before it can be cropped
library(lwgeom)
st_is_valid(map_1300)
map_1300_fixed <- st_make_valid(map_1300)
st_is_valid(map_1300_fixed) # this doesn't work
sf::sf_use_s2(FALSE) # this does work
map_1300_crop <- st_crop(map_1300, bbox_cities_sf)
sf::sf_use_s2(TRUE) # putting settings back to their default
# plot the cropped borders and cities
plot(map_1300_crop[12], reset=F)
plot(cities_165_crop[4], add=T)
# Make the figure for the paper...
state_map <- tm_shape(map_1300_crop) + tm_borders(col = "black")
state_map
mycols <- colors()[c(142, 54, 35, 41, 24)]
cities_map_color <- tm_shape(cities_165_crop) +
tm_bubbles("mortality", col = "mortality", border.col= "black",
border.alpha = 0.5, style="fixed",
breaks = c(-Inf, 20, 40, 60, 80, Inf),
palette=mycols, title.size = "Mortality Rate (Size)",
scale = 2.0, title.col ="Mortality Rate (Color)") +
tm_layout(legend.stack = "vertical", legend.text.size = .80,
frame = F, inner.margins = .10,
legend.position = c(-.15, 0.05)) +
state_map
cities_map_color
# save map
tmap_save(cities_map_color, filename="cities_map_color.png",
height=8.5, width=11, units="in", dpi=300)
# Code to create Figure A1 in "Pandemics and Cities: Evidence from the Black Death and the Long-Run"
# Author: Noel Johnson
# Date Created: 5-14-20
# Last updated: 9-5-23
library(sf)
library(raster)
library(tidyverse)
library(tmap)
library(haven)
# Need the following features in the map:
# (1) Modern Country Borders
# (2) Cities: ASTRAHAN, FEODOSIJA, MESSINA, GENOVA, IZMAIL, WIEN, MOSKVA, NOVGOROD, VISBY, Prague, Leipzig, Halych
# (3) Rivers: Danube, Dniester, the river from Astrakhan to Moscow, moscow to Novgord
# (4) Medieval Trade Routes: the medieval land route (Via Regia),
# (5) seas
# coordinates of the cities not in Bosker (from wikipedia)
# c("ASTRAHAN",51.166667, 71.433333)
# c("FEODOSIJA",45.048889, 35.379167)
# c("VINICA",45.351667, 28.836389)
# c("MOSKVA",55.755833, 37.617222)
# c("NOVGOROD",58.55, 31.266667)
# c("VISBY",57.634722, 18.299167)
# set working directory
setwd("/Users/noeljohnson_laptop/Dropbox/Research/Plague City Growth/JUE RR1 Nov 2022/figures_replication/Figure_A1")
# bring in the cities
cities_raw <- read_dta("All_Cities.dta")
# keep the variables you want
cities_raw <- cities_raw %>% select(city_jjk, longitude, latitude)
# keep the cities you want
target <- c("WIEN", "PRAHA", "LEIPZIG")
cities_raw <- cities_raw %>% filter(city_jjk %in% target)
# add the cities not contained in Bosker that you acquired coordinates for form Bosker
cities_raw <- cities_raw %>%
rbind(c("ASTRAHAN (1345)",48.05,46.35)) %>%
rbind(c("KAFFA (1346)",35.379167,45.048889)) %>%
rbind(c("VINICA",28.836389,45.351667)) %>%
rbind(c("MOSKVA",37.617222,55.755833)) %>%
rbind(c("NOVGOROD",31.266667,58.55)) %>%
rbind(c("VISBY",18.299167,57.634722)) %>%
rbind(c("HALYCH",24.728611,49.124722)) %>%
rbind(c("GENOVA",8.946256,44.4056499)) %>%
rbind(c("MESSINA (1347)",15.5540152,38.1938137))
# clean up the cities data
glimpse(cities_raw)
cities_raw$longitude <- as.numeric(as.character(cities_raw$longitude))
cities_raw$latitude <- as.numeric(as.character(cities_raw$latitude))
glimpse(cities_raw)
# spatialize the cities data
cities <- st_as_sf(cities_raw, coords=c("longitude", "latitude"), crs=4326)
# take a look at the data
cities
glimpse(cities)
head(cities)
plot(cities[2], col="red", pch=16, reset=F)
# Make a bounding box of the cities so you can create one a little bigger
bbox_cities <- st_bbox(cities, crs=4326)
bbox_cities
xrange <- bbox_cities$xmax - bbox_cities$xmin # range of x values
yrange <- bbox_cities$ymax - bbox_cities$ymin # range of y values
bbox_cities[1] <- bbox_cities[1] - (0.20 * xrange) # xmin - left
bbox_cities[3] <- bbox_cities[3] + (0.20 * xrange) # xmax - right
bbox_cities[2] <- bbox_cities[2] - (0.20 * yrange) # ymin - bottom
bbox_cities[4] <- bbox_cities[4] + (0.20 * yrange) # ymax - top
# crop the cities using your bounding box
cities <- cities %>%
st_crop(bbox_cities)
# Bring in medieval roads shape file
med_roads <- st_read("Europe_Medieval_Trade_Routes/MedRdsProj.shp") %>%
st_transform(4326)
# take a look at the data
glimpse(med_roads)
head(med_roads)
# get rid of everything except the gis field
med_roads <- med_roads %>% dplyr::select()
# crop the roads
med_roads <- med_roads %>%
st_crop(bbox_cities)
glimpse(med_roads)
# add Via Reggia
crs = st_crs(4326)
s1 <- rbind(c(24.728611,49.124722),c(14.4378005, 50.0755381))
ls <- st_linestring(s1)
b = st_sf(geometry = st_sfc(ls), crs = crs)
med_roads <- med_roads %>%
rbind(b)
# plot the roads
plot(med_roads[1], col="brown", reset=F, add=T)
# Bring in rivers shape file: round 1
rivers <- st_read("Rivers/MajorRivers.shp") %>%
st_transform(4326)
glimpse(rivers)
rivers <- rivers %>% dplyr::select(NAME)
rivers <- rivers %>%
st_crop(bbox_cities)
glimpse(rivers)
plot(rivers[2], col="blue", reset=F, add=T)
# Bring in rivers shape file: round 2
rivers2 <- st_read("Rivers/Rivers_conic.shp") %>%
st_transform(4326)
glimpse(rivers2)
rivers2 <- rivers2 %>% dplyr::select(name)
rivers2 <- rivers2 %>%
st_crop(bbox_cities)
glimpse(rivers2)
target <- c("Dniester", "Volga", "Oka")
rivers2 <- rivers2 %>% filter(name %in% target)
plot(rivers2[2], col="blue", reset=F, add=T)
# Bring in the modern country borders
modern_countries <- st_read("ModernCountries/Longitude_Graticules_and_World_Countries_Boundaries-shp/99bfd9e7-bb42-4728-87b5-07f8c8ac631c2020328-1-1vef4ev.lu5nk.shp") %>%
st_transform(4326)
# Clip the countries using bounding box based on the cities data
modern_countries <- modern_countries %>%
st_crop(bbox_cities)
glimpse(modern_countries)
head(modern_countries)
modern_countries <- modern_countries %>% dplyr::select(CNTRY_NAME)
glimpse(modern_countries)
plot(modern_countries[2], col="NA", reset=F, add=T)
# Bring in seas shape file
seas <- st_read("Seas/ne_10m_ocean.shp") %>%
st_transform(4326)
glimpse(seas)
# there is a ring intersection error in the seas .shp file. So need to fix the geometry.
sf::sf_use_s2(FALSE) # this does work
seas_crop <- st_crop(seas, bbox_cities)
sf::sf_use_s2(TRUE) # putting settings back to their default
# make figure A1
state_map <- tm_shape(modern_countries) + tm_borders(col = "grey")
rivers_map <- tm_shape(rivers) +
tm_lines(col="blue", lwd=3) +
tm_add_legend("line", col = "blue", title = "Major Rivers", lwd=3)
rivers_map2 <- tm_shape(rivers2) +
tm_lines(col="blue", lwd=3)
med_map <- tm_shape(med_roads) +
tm_lines(col="black") +
tm_add_legend("line", col = "black", title = "Trade Routes")
seas_map <- tm_shape(seas_crop) +
tm_polygons(col="grey")
cities_map <- tm_shape(cities) +
tm_text("city_jjk", size=0.75,ymod=1,xmod=.75, fontface="bold") +
tm_dots(col="red",size=.40) +
tm_add_legend("symbol", col = "red", title = "Trade Cities")
full_map <- rivers_map + rivers_map2 + med_map + cities_map + seas_map
full_map <- full_map + tm_layout(legend.position = c(.77,.70),legend.frame=T)
full_map
# save the map (note: the orange arrows were inputted by hand)
tmap_save(full_map, filename="Trade_Figure1.png",
height=8.5, width=11, units="in", dpi=300)
# Code to create Figure A10 in "Pandemics and Cities: Evidence from the Black
# Death and the Long-Run"
# Author: Noel Johnson
# Date Created: 10-17-19
# Last updated: 9-8-23
library(raster)
library(sf)
library(tidyverse)
library(tmap)
setwd("/Users/noeljohnson_laptop/Dropbox/Research/Plague City Growth/JUE RR1 Nov 2022/figures_replication/Figure_A10")
# Load in the data sets, make spatial, and project
# Define the projection
EEC <- "+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=43 +lat_2=62 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs"
# Load in the cities data
cities_full <- read.csv("list_467_165_cities.csv")
head(cities_full)
# spatialize the cities data
cities_full_spatial <- st_as_sf(cities_full,
coords = c("longitude", "latitude"), crs = 4326)
# project the cities data
cities_165_proj <- st_transform(cities_full_spatial, EEC)
st_geometry(cities_165_proj)
# Load in the modern country borders
modern_borders <- st_read("ModernCountries/Modern Europe Projected.shp")
# project the modern country borders
modern_borders_proj <- st_transform(modern_borders, EEC)
st_geometry(modern_borders_proj)
# need to crop the maps
# original Extent of modern countries
# c(xmin=-2639987, xmax=1655719, ymin=3243609, ymax=8097861))
modern_borders_clip <- st_crop(modern_borders_proj,
c(xmin=-1100000, xmax=1655719,
ymin=3243609, ymax=6800000))
# Make a simple plot of the cities and borders to make sure it looks right
plot(modern_borders_clip["name"], col=NA, reset = FALSE)
plot(cities_165_proj[5], add = TRUE)
# Create two data sets:
# (1) 467 cities with pop >1000 in 1300;
# (2) 165 cities in our sample (mort = 1)
head(cities_165_proj)
pop1300_cities <- cities_165_proj %>%
select(-samplemort)
sample_cities <- cities_165_proj %>%
filter(samplemort==1) %>%
select(-samplemort)
# Make the figure for the paper
state_map <- tm_shape(modern_borders_clip) + tm_borders(col = "black")
breaks <- c(50, 100, 150, 200, 250)
cities_map <- tm_shape(pop1300_cities) +
tm_bubbles("pop1300", "black", border.col = "white",
shape = 1,
sizes.legend = breaks,
title.size ="Population in 1300",
scale = 1.75, perceptual = F, sizes.legend.labels = c("50", "100", "150", "200", "250")) +
tm_layout(legend.stack = "vertical", legend.text.size = .75, frame = F, inner.margins = .10, legend.position = c(-.05, 0.36)) +
state_map
full_map_x <- tm_shape(sample_cities) +
tm_dots("black", size = .4, shape = 4) +
tm_add_legend(col = "black", type = "symbol", shape = 4, title = "Sample Cities") +
cities_map
full_map_x
# save the map
tmap_save(full_map_x, filename="figure_a10.png",
height=8.5, width=11, units="in", dpi=300)
